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The nature of electron correlations in bilayer graphene has been investigated. An analytic expres- 
sion for the radial distribution function is derived for an ideal electron gas and the corresponding 
static structure factor is evaluated. We also estimate the interaction energy of this system. In 
particular, the functional form of the pair-correlation function was found to be almost insensitive to 
the electron density in the experimentally accessible range. The inter-layer bias potential also has a 
negligible effect on the pair-correlation function. Our results offer valuable insights into the general 
behavior of the correlated systems and serve as an essential starting-point for investigation of the 
fully-interacting system. 
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Despite intense studies for many decades, the im- 
portant role of many-particle correlations in electron 
liquids 1 , particularly in systems with reduced dimen- 
sions, remains a challenging issue in condensed-matter 
physics. This subject has become even more pressing 
in recent years as the physical properties of graphene 
have been unmasked at a rapid pace^. Monolayer and bi- 
layer graphene are totally new classes of two-dimensional 
electron systems with unusual band structures and chiral 
charge carriers. The influence of electron correlations on 
various physical properties of the chiral two-dimensional 
electron gas in monolayer graphene has been the subject 
of several investigations 3-5 . These were all carried out 
with the aid of established pertubative techniques. On 
the other hand, our earlier work based on exact analyt- 
ical treatment^ indicated that the electron correlations 
completely vanish from the two-particle kinetic energy of 
monolayer graphene, a fact which was attributed to the 
specific spinor structure of the single-particle wave func- 
tions which in turn is a direct manifestation of the chiral- 
ity of the massless Dirac fermions in monolayer graphene. 
Such cancellations were not found to occur for electrons 
in bilayer graphene^, due to the massive chiral nature of 
the low-energy quasiparticles. Clearly, we have a long 
way to go in order to find a satisfactory understanding 
of the role interactions play in these unique electron sys- 
tems, but it is evident that the effects of correlations in 
bilayer graphene is an important and relevant issue. 

In this Rapid Communication, we will lay out the 
foundation for the process of establishing the behavior 
of the correlation function of interacting electrons in bi- 
layer graphene, which is an essential step for evaluation of 
the thermodynamic properties of this system. We derive 
an analytic expression for the pair-correlation function 
(PCF) of an ideal electron system, and use it to compute 
the corresponding static structure factor as a function 
of the electron density. We make a detailed compari- 
son of the PCF with the same quantity in a traditional 
two-dimensional electron gas (2DEG), and compute the 
exchange energy for the bilayer graphene system. Evalu- 
ation of the PCF with full electron correlations included 
is certainly a very arduous task£, and has not yet been 



attempted for bilayer graphene. Our present approach 
is an important and necessary first step in describing a 
fully correlated system, but it already provides valuable 
insights into the general behavior of these functions. 

Bilayer graphene has a hexagonal Brillouin zone where 
the low-energy features are located near two inequivalcnt 
corners, called the K points^. There are four low-energy 
7r bands for each spin species in the vicinity of each K 
point. In unbiased bilayer graphene the second and third 
of these bands (called the 'low-energy branches') touch 
exactly at the charge-neutrality point, making intrinsic 
bilayer graphene a zero-gap semiconductor. There are 
two more bands (the 'split' branches), each seperated 
from the low-energy branch by the inter-layer coupling 
parameter j 1 ss 0.4 eV. We label these bands as follows: 
The conduction and valence bands are given by v = c, v; 
the branches are b — l,s; the valleys are £ = K, K'\ and 
the spins are a =j, |. Adding the electron wave vector 
k, we have the complete set of quantum numbers A = 
{za , b x , k A , £ A , <t a }. At half-filling, all eight valence bands 
are filled, and the eight conduction bands are empty. 

Spin- and valley- dependent contributions to the PCF 
are defined asi 
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where ^ „ x ^ x ( r ) = Sk A v b x ^A a A 1S ^ ne ne ^ operator for 
an electron in valley £a with spin er A . The total PCF can 
be expressed in terms of these functions as (7(1*1, ra) = 
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To evaluate this expression we substitute the well- 
known form of the single-particle wave functions in bi- 
layer graphene 
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where s A and S A are respectively the valley and spin 
parts of the wave function, A is the normalization area, 
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and the functional form of the wave function components 
and single-particle energies are easily derived from the 
Schrodinger equation for the tight-binding Hamiltonian. 
Substituting Eq. © into Eq. flU, we find that 
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carriers (also called the excess density) n cc = n — n 
which may be either positive (for electrons) or negative 
(for holes). Then, the sums over occupied wave vector 
states must be taken independently for each combina- 
tion of band and branch quantum numbers. Taking the 
limit of an infinite system (with the electron density held 
constant) means that we can replace the sums over wave 
vectors with two-dimensional integrals. The integrals 
which result from this procedure are not automatically 
convergent for large wave vectors. Therefore, we must 
introduce a cut-off wave vector using some physical rea- 
soning. Consideration of the lattice structure shows that 
each unit cell contributes four tt electrons (one per car- 
bon atom), so that the density of electrons at half- filling 
is n a = 8/ (VSa 2 ) where a s» 2.46A is the lattice constant. 
Therefore, we can set the wave vector cut-off A because 
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with r = r x — r 2 and where / A is the occupancy of state 
A. As expected, the off-diagonal components of the PCF 
are constant with unit value. We have also assumed 
that n o- A £ A ( r i) = n /4 (i- e - * na t electrons are equally 
distributed between the valley and spin components and 
that electron density is uniform in space). 

To procede, we must be careful about how we define 
the various densities. The total density of electrons is 
denoted by n, but we also consider the density of charge 

I 

As an example, in intrinsic graphene (where the Fermi energy is exactly at the charge neutrality point) the valence 
bands are all filled and the conduction bands are all empty. Therefore, the sum over bands and branches in Eq. ([3]) 
becomes 
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where Nq is the total number of electrons at half-filling. 
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Using the expressions for the wave functions in Eq. ([2]) and evaluating the elementary integrations over the angles, 
we arrive at 
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where all terms are evaluated with v = v, and Jo(^) and 
Ji(x) are the zeroeth order and first order cylindrical 
Bessel function respectively. In the case of positively- 
doped graphene (where the charge carriers are holes and 
the Fermi energy is in the valence band), we assume 
that for moderate densities only the low-energy band 
is depopulated^ so that the lower integration limit be- 



comes the Fermi wave vector k F = ^nn cc when b = I. 
For negatively-doped graphene, each squared term in 
Eq. Q gains a contribution from the low-energy con- 
duction band (y = c, b — I) with the Fermi wave vector 
replacing A as the upper limit in this integral. 

In order to obtain the PCF, these integrals are evalu- 
ated numerically and the resulting function is plotted in 
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FIG. 1. a) The pair-correlation function for several values 
of the electron density density. Black solid line: n cc = 0; 
red dashed line: n cc = 10 13 cm~ 2 ; green dash-dot line: n cc = 



5 x 10 cm ; blue dotted line: 



= lO^crrT 



b) Pair- 



correlation function near the exclusion hole edge. Densities 
are as in a), c) Comparison with 2DEG PCF. Red dashed 
line: 2DEG n = 10 13 cm~ 2 ; black solid line: 2DEG n = 7.6 x 
10 15 cm~ 2 ; black dotted line: bilayer n cc = 0. 



Fig. Q] for various densities. The behavior of the function 
is clearly similar to that in a conventional 2DEG^i, with 
an exchange hole with radius approximately 5A. The rea- 
son that g(0) is finite is explained as follows: The PCF 
evaluated here specifies all-but-two quantum numbers. 
Therefore for any given combination of band and branch, 
there can be an electron at r = with one of three other 
combinations which does not violate the Pauli principle. 
Hence the minimum value of the PCF is 3/4, as seen in 
Figffl If we were to calculate the g(r) for a fully-specified 
combination of valley, spin, band and branch then this 
function would indeed go to zero at the origin, just as it 
does for the conventional 2DEG. 

The dependence of the PCF on the density is tiny 
for physically reasonable values of the excess density. 
The reason for this tiny variation is that the electrons 
in the filled valence bands contribute more to the sum 
over states than those in the partially filled conduction 
band. The PCF contains essentially an average over all 
particles -k ^ (where N is the total number of elec- 
trons and i runs over all filled states). The intrinsic 
density of electrons due to the valence bands which are 
filled in the charge- neutral case is no ~ 7.6 x 10 15 cm~ 2 , 
which is much greater than the density of charge carri- 
ers n cc < 10 14 cm -2 due to the excess density induced 
by gating or doping. Therefore when the average over all 
states is taken, the effect of the partially filled conduction 
band (or partially empty valence band) is swamped by 
the contribution from the filled valence band. This effect 
is highlighted by comparison with the non-interacting 
PCF in a traditional semiconductor 2DEG (in the lower 
inset to Fig. [J). When the 2DEG PCF is plotted for 
n = 1 x 10 13 cm -2 , the exchange hole is much larger than 
in graphene. But when the total density 7.6 x 10 15 cm -2 
is used the exchange hole is of a much more similar size. 

Once we obtain the radial distribution function, the 



FIG. 2. The static structure function for the same densities 
as in Fig. QJ,). The units are r = hvp/^i « 1.65nm. Main 
plot: The whole wave vector range. Upper-left inset: The 
low- wave vector region. Lower-right inset: The wave vector 
region near the cut-off. 



static structure factor for the system can be derived from 
the following expression 1 
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which is, in principle, an experimentally observable func- 
tion via X-ray and neutron diffraction where the corre- 
lation functions are usually extracted from the measured 
diffraction intensity profile 12 . We have evaluated the in- 
tegral numerically and the resulting function is plotted 
in Fig. [2] for several values of the electron density. We 
see that the variation with density is rather small, but 
at low wave vector, S(k) increases with density (upper- 
left inset) while at high wave vector the opposite is true 
(lower-right inset). The structure factor is almost lin- 
ear even up to the wave vector cut-off A. This behav- 
ior has been noticed before in the context of monolayer 
graphene^i This is noticably different from the result 
for the conventional 2DEG 1 *, where the static structure 
function is roughly linear at small wave vector but satu- 
rates at 5 I 2DEg(' c ) = 2 f° r ^ > 2fc F . We emphasize that 
the static structure function of bilayer graphene behaves 
similarly to the conventional 2DEG at small wave vector, 
but like monolayer graphene at large wave vector. This 
behavior might be expected when the quadratic-to-linear 
crossover in the hyperbolic band structure is considered. 

Finally, we calculate the exchange energy per electron 
associated with the exchange-correlation hole 
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where V(r) is the Coulomb potential and we use the full 
g(r). This function is linear in the quasi-particle density, 
with £7(0) -2.5 eV. 

Let us now turn our attention to the effect of a fi- 
nite inter-layer bias potential on the radial distribution 
function. When an electrostatic potential is applied per- 
pendicularly to the plane of the graphene, a gap opens 
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FIG. 3. The change in g(r) with the inter-layer bias potential 
for several densities. Main plot: U = 500 me V; inset: U = 
100 meV. All densities are in cm~ 2 and line styles are the 
same as in Fig. [ij,). 

at the charge-neutrality point, and the shape of the low- 
energy bands changes to a 'Mexican hat' form^. This 
also changes the form of the wave functions, and causes 
the Fermi surface to become ring-shaped for small charge 
carrier density^. Therefore, the integration limits in 
Eq. (01 change if the Fermi energy E F < U/2. In that 
case, integrals relating to partially-filled bands become 
f* F dk — !• dk with 
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On the other hand, if E F > U/2 [which occurs when 
k F > U/(hv F )] then k F = ^/im cc as before. 

We plot the change in the PCF with the introduction 
of a bias Ag = gjj(r) — 5o( r ) as a function of the inter- 



particle separation in Fig. [3] We see that the change 
is greatest at small charge carrier density and large U . 
However, overall the change is very small which is pre- 
dictable since the PCF is related to the electron wave 
functions, and the inter-layer potential only induces a 
change for E ~ U <C E(A). Similarly, the static structure 
factor shows only very small deviation from the U = 
results for finite U . 

In conclusion, we have investigated the PCF and the 
corresponding static structure function for an ideal gas of 
electrons in bilayer graphene and compared it to the same 
quantity in the traditional 2DEG system. We have found 
behavior quite similar to that of the conventional 2DEG 
at equivalent density, in that an exchange hole is formed 
with density-dependent radius. However, the manifesta- 
tion of effects due to the bands, especially the existence 
of the filled valence band means that the dependence of 
these functions on the density of charge carriers is min- 
imal in the experimentally accessible range. We have 
evaluated these functions for the gapped system as well, 
and found that the effect of the inter-layer bias poten- 
tial on these quantities was also negligable. This general 
picture will also be true for all Dirac-like systems which 
have filled valence bands. In the case when the many- 
body correlations are taken into account, we expect very 
similar behavior for the dependence on density because 
electron-electron interactions do not alter the situation 
of filled valence bands. In monolayer graphene, we pre- 
viously expected that the functional form of the PCF 
remain insensitive to electron density in order to explain 
the observed behavior of the electron compressibility^. It 
is interesting to observe a similar situation in the case of 
another Dirac-like graphene system based on very general 
considerations. 
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